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Treating dynamical stability as an observable: a 5:2 MMR 
configuration for the extrasolar system HD 181433 
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ABSTRACT 

The three-planet extrasolar system of HD 181433 has been detected with HARPS. 
The best-fit solution, announced by the discovery team, describes a highly unstable, 
self-disrupting configuration. In fact, a narrow observational window, only partially 
covering the longest orbital period, can lead to solutions representing unrealistic sce- 
narios. Taking into account the dynamical stability as an additional observable while 
interpreting the RV data, we can analyse the phase space in a neighbourhood of the 
statistically best-fit and derive dynamically stable configurations that reproduce the 
observed RV signal. Our Newtonian stable best-fit model is capable of surviving for 
at least 250 Myrs. The two giant companions are found to be locked in the 5:2 MMR 
as Jupiter and Saturn in the Solar System. This mechanism does not allow close en- 
counters even in case of highly eccentric orbits. Moreover, planets c and d are located 
in regions spanned by many other strong low-order MMRs. We study the dynamics 
of some plausible scenarios and we illustrate the behaviours caused by secular apsidal 
resonances and mean motion resonances. Furthermore, we find a terrestrial planet in 
the habitable zone of HD 181433 can retain stability. Apart from filling an empty gap 
in the system, this body could offer a harbour for life indeed. Additional measurements 
are necessary in order to investigate this hypothesis and can confirm the predictions 
outlined in the paper. 

Key words: planetary systems, dynamical evolution and stability, celestial mechan- 
ics, radial velocity technique - stars: individual: HD 181433 - methods: N-body sim- 
ulations, numerical, statistical. 



1 INTRODUCTION 

Nowadays, more than 50 multi-planet systems are known to 
harbour 2-6 planet^]. The recent years have seen a prolif- 
eration of multiple-planet systems thanks to the increase in 
both instrument precisions and duration of several planet 
search programs. This has allowed the detection of longer 
periods planetary signatures, as well as planetary signatures 
with lower amplitude. 

An increasing sample of found planets improves our 
knowledge of their distribution in the mass-period diagram 
and allows comparison with theoretical predictions (e.g. 
Wright et al. 2009). In fact, dynamical analysis of planetary 
systems can both precisely determine their orbital architec- 
tures and constrain their evolutionary histories providing a 
test bed for planetary formation and evolution theories. To 
make these investigations we need measured orbital param- 
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eters as accurate as possible. Unfortunately, the accuracy 
of such measurements are limited due to uncertainties and 
degeneracies inherent to the Radial Velocity (RV) discovery 
technique. Nevertheless, the RV method is the most efficient 
technique for detecting extrasolar planets with more than 
90% of all currently known planets being either detected or 
characterized using this method. 

As the time baseline becomes larger, it is possible to 
distinguish a trend in the RV signals due to long-period 
outer companions that have not completed a single orbit 
yet. In this case, it can happen either that the profile of 
X 2 is very smooth, or that it does not have a well definite 
minimum, as a result the confidence levels may comprise 
large intervals of the parameters fitted. 

The RV signal does not provide any direct information 
on the real masses and the orbital inclinations of the planets; 
even if we make the assumption the system is coplanar and 
seen edge-on, an N-planet configuration is described by some 
5N+1 parameters. The a priori unspecified number of plan- 
ets, narrow observational windows, stellar jitter and weakly 
constrained orbital parameters, can lead to not unique so- 
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lutions and (quite often) to best-fits representing unreal- 
istic scenarios. In fact, these solutions can present well- 
constrained minimum masses but also a poorly constrained 
eccentricity for the outermost planet; i. e., in the statistically 
optimal best-fit solutions, the eccentricities can be large and 
can rapidly (on the timescale of thousa nds of years) bring 



to ca tastrophic collisional instabilities (jGozdziewski et al 
l2008h . 

According to the Copernican assumption that we are 
not observing the universe at a privileged time, the detec- 
tion of a rapidly unstable system during a few years of RV 
observations is not likely, then a fit which corresponds to a 
quickly unstable configuration is also doubtful. Therefore, 
we can aim to put limits on the masses and orbital ele- 
ments of the planets by investigating and finding the plau- 
sible and stable solutions. In this logic, the dynamical sta- 
bility is an additional observable that must be taken into 
consideration when interpreting the RV data. It turns into 
a discriminating element especiall y when the longest orbital 
period is only partially covered ()Murrav fc Holmanl 1200 ll : 
IGozdziewski et alJl2003 ). 

Often, a stability criterion is applied once the best-fit 
unstable solution is found, subsequently the orbital elements 
are tuned to get a stable configuration. However this ap- 
proach does not necessarily provide stable fits that are si- 
multaneously optimal in term of \ 2 or RMS. Most of the 
times, with the term stable we indicate a configuration which 
does not disrupt or change qualitatively during a period of 
time of the order of million years. The literature is plenty of 
studies which take into account stability when modelling the 
RV da ta. Here we just mention the work of: IVeras fc Fordl 
(2010) reporting on the system stability, secular evolu- 
tion and the extent of the resonant i nteractions for 5 dy - 
namically active multi-planet systems; Wrig ht et al.l (|2009I ) 
which derive updated orbital parameters for a number of 
systems considering mutual in t eract ions between planets; 
IGozdziewski et al.1 (|2003l . 120061 . |200S| ) which directly elim- 
inate unstab le solutions during t he fitting procedure using 
GAMP; and ICorreia et al.1 (|2010h which give constrains on 
the inclination with respect to the line of sight for some of 
the planets in the GJ 876 system. We also should remind the 
study about the directly imaged system HR 8799, where the 
difficulties in finding regions of stable motion may indicate 
the system is undergoing a phase o f planet-planet scattering 
IGozdziewski fc Migaszewskill2009l ). 

Thanks to powerful instrumentations like the HARPS 
spectrograph at La Silla Observatory, the Radial Veloc- 
ity acc uracy is increasin g rapidly breaking the barrier of 
1 m/s (|Pepe et al.l 120031 ') . This has permitted the detec- 
tion of weaker signals consenting the discovery of some 
of the lowest-mass planets identi fied such as: GJ 8 76 d 
(iRiveraet al.ll2005l). H P 40307 b (iMavor et al.1 1200981). 6 1 
Vir b (|Vogt et al.ll2010l ) and GJ 581 e l|Mavor et al.ll2009bl) . 

The planetary system of HD 181433 has been discov- 
ered with HARPS. It has been reported to contain 3 plan- 
ets: two Jupiter-clas s planets and a Super-Earth of 7.5 m© 
l|Bouchv et al.ll2009r i. 

Inspired by the peculiar properties of the system, which 
includes two giant planets and one rocky planet all in high 
eccentric orbits, we aimed to study the past and future evo- 
lution of the system. Unfortunately, the published configu- 



ration is unstabkfl The model in which the initial eccentric- 
ities of the planets are reduced by one sigma quickly leads to 
disruption too. The fate does not change when we assume, in 
addition, a mutual orbital inclination of 20 ° between planets 
c and d. 

These attempts evidence the necessity of doing an anal- 
ysis from scratch in order to get a self-consistent solution 
compatible with the data. In this paper we examine the 
availa ble RV data of the HD 181433 system (|Bouchv et al.l 
2009) taking a more general approach, going beyond a for- 
mal fit of the Keplerian orbital elements. Even if the RV 
observations do not span a single period of the outermost 
planet, we can give reasonable constraints on the orbital el- 
ements of the poorly sampled third planet by studying the 
dynamics of the system. In the Keplerian fit this important 
information is completely omitted. 

In many situations, the interactions between planets are 
negligible and can be ignored. So the RV signal is just a lin- 
ear superposition of different Keplerian RV curves. On the 
other hand, planetary interactions and resonances can be im- 
portant and must be taken into account performing numer- 
ical integrations. In the short time scale, these interactions 
can cause significant variations in the orbital parameters of 
the planets. As in our case, an ensemble of constant Keple- 
rian orbital elements is not adequate to model the RV data 
and an A^-body Newtonian model should be applied. 

However, it can also happen that good Newtonian fits 
to the data produce planetary orbital parameters that are 
stable for the period of observations but lead to disruption 
on timescales substantially shorter than the age of the plan- 
etary system. In this respect, long-term stability is an ad- 
ditional but necessary constraint that must be satisfied by 
multiplanet fits. 

In Section [2] we perform an independent analysis of the 
RV data for HD 181433. We probe the phase space of the or- 
bital parameters looking for likely configurations stable for 
long timescales, say millions of years. We assume the motion 
is described by Newtonian interacting orbits. In Section [31 
we present the dynamical study of the stable best-fit solu- 
tion and we analyse the behaviour of other plausible stable 
configurations. In particular, we focus on the description 
of secular apsidal resonances (SARs) and mean motion res- 
onances (MMRs). In Section 3] we briefly summarize our 
findings, we discuss on the possibility of a terrestrial planet 
in the habitable zone and we make some predictions about 
what we may expect from further observations. 



2 RADIAL VELOCITY DATA ANALYSIS 



iBouchv et ail (2009) announced the detection of 3 planets 
around HD 181433, a K3IV star, considering 107 RV mea- 
surements which covered more than 4 years, from June 2003 
to March 2008. The median uncertainty for the RV is 0.53 
m/s with most of the data in the range 0.4-1.0 m/s. The 
peak-to-peak velocity variation is 48.12 m/s, while the ve- 
locity scatter around the mean RV in the measurements is 
13.86 m/s. 

The data do not completely cover a full period for the 
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third planet. In fact, what is possible to spot i t is just an ad- 
dition al long-term trend which is modeled bv lBouchv et al.l 
(2009) as being produced by a planet of minimum mass 
msini = 0.54 mj up on a Keplerian orbit with a period of 
about 6 years and e = 0.48. This model is unstable on the 
order of just thousands of years. The instability arises due 
to close encounters between planet c and planet d. 

We perform a r e-analysis of the HAR PS data using the 
Systemic Consol^fl l|Meschiari et al.ll2009l ). Systemic has al- 
re ady been u s ed to deri ve orbital fits in other works such 
as lVogt et all (|2010h and lMeschiari et alj |201ll ). The list of 
available tools include: Lomb-Scargle (LS) periodogram to 
identify periodicities in the RV dataset, Lomb-Scargle pe- 
riodogram of residual to study periodicities in the residual 
RV dataset, simulated annealing for global multiparameter 
optimization, while for local multidimensional optimization 
there are the Levenberg-Marquardt (L-M) scheme which en- 
sures a rapid convergence an d the Nelder-Mead (sometimes 
called AMOEBA) algorithm (jPress et al.lll992T ). 

We held th e stellar mass fix ed, adopting the value 
M» = O.781M l|Sousa et al.ll2008f ). We make the assump- 
tion the system is coplanar and viewed edge-on. This con- 
jecture diminishes the quantity of potential orbital configu- 
rations greatly, but the plane (ad, e.d) is dynamically repre- 
sentative for_th£_sj^stejnJ^thesense that it crosses all res- 
onances (|Robutel fc Laskarlliobll ). When long enough time- 
series of precision data are available, the effects of mutual 
interactions part of the Newtonian model can potentially 
help in determining or estimating the masses and inclina- 
tions for the planets (see the discussion about this case in 
Section [2~2)l. 



2.1 The Keplerian 3-planet best-fit 

The LS periodogram shows a peak at P — 1171.35 days with 
an estimated false alarm probability (PAP) of ~ 2 x 10 -19 . 
Fig. Q] shows the LS periodogram of the full RV data set and 
the periodogram of sampling. The latest shows peaks that 
are related with periodicities in the cadence of observations, 
for instance these can arise from the solar and sidereal day, 
the synodic month and the solar year. 

The residuals periodogram reveals an additional signal 
at P = 9.37 days with FAP of « 3.8 x 10 -5 . The best 2- 
planet Keplerian fit yields residuals with an rms scatter of 
2.44 m/s and reduced chi squared Xred = 15.7. The jitter 
for HD 181433 i.e. the jitter required to have the Xred equal 
to 1.0, is 2.35 m/s. Fig. illustrates the periodogram of 
residuals to the two-planet solution. 

To model this long-term trend, we make the starting 
guess of a planet in an outer 2:1 resonance with planet c, we 
adjust the mass to match the amplitude of the signal and set 
a small eccentricity. At this point, a Keplerian fit using the 
L-M algorithm naturally evolves to a solution compatible 
with the one bv lBouchv et alJ l|2009l ). The best 3-planet fit 
achieves a Xred = 4.6 with an rms scatter of 1.34 m/s and 
expected jitter of 1.17 m/s. We are not aware of how Bouchy 
et al. have obtained a lower Xred an d a lower value for rms 
for the same solution. 



3 Available at http : / /www . oklo . org 
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Figure 2. Residuals periodogram to the two-planet Keplerian fit. 

The left panel of fig.[3]shows the best- fit orbital configu- 
ration at the epoch of the first observation BJD 2452797.87. 
The orbits of planets c and d cross each other, because of 
the strong interactions collisions/ejections occur. 

2.2 The Newtonian 3-planet stable best-fit 

Following the argument of lAnglada-Escude et all (|2010h 
that eccentric orbital solutions can mimic the signal of two 
planets in 2:1 resonant orbits, we have also tested the hy- 
pothesis of a planet in an inner 2:1 resonance with planet c 
but it was not possible to achieve any significant improve- 
ment to the goodness of the fit with respect to the 2-planet 
solution. 

The problem of exploring a 16-parameter phase space 
with stability as additional requirement, can get a first sim- 
plification by arguing that the elements of the inner plan- 
ets are well constrained by observations. In fact, even if we 
set different starting points for their parameters, the fits 
for them converge substantially to the same values as these 
signals are well sampled. A confirmation to this argument 
comes from the eccentricity of planet c, e c , which is a very 
discriminating parameter toward the stability of the system: 
if we try to constrain e c to lower values the fitting we achieve 
is poor. 

Concerning the parameters that describe planet d, we 
notice higher values for the eccentricity are preferred by the 
fitting. Therefore, at the end the problem can be reduced 
in finding for each reasonable Pd the largest value for for 
which planets c and d do not undergo instability. Likewise, 
we can say once we have a stable solution for a pair (Pd 
- ed) we want to investigate if it possible to get a differ- 
ent pair which generates a stable configuration having the 
same or lower Xred- We can describe that as being an empir- 
ical Bayesian approach of inferring the stable best-fit rather 
than a frequentist approach which involves a time consuming 
number of simulations. Here the investigation is conducted 
by evidences like the ones given by the collision line (see 
later on in this Section) and the outcomes of previous sim- 
ulations. 

The second and third planet reside in regions spanned 
by a number of strong low-order MMRs (see Fig. [5] later 
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Figure 1. Periodograms for HD 181433. Left: LS periodogram of the RV data, a peak at P = 1171.35 days is present. The three 
horizontal lines represent, from top to bottom, the 0.1 %, 1.0 % and 10.0 % FAP levels, respectively. Right: Periodogram of sampling 
showing peaks that are related with periodicities in the cadence of observations. 




Figure 3. Orbital views for the HD 181433 system, the position of the planets along their orbits is the one at the first epoch of 
observation. The straight lines point toward the pericentres. The osculating elements are valid for the first epoch of observation. Left 
Panel: Best fit solution, the orbits of planets c and d cross each other and collisions/ejections occur. Right Panel: Stable best fit solution 
planets c and d are in 5:2 MMR. 



on). We are aware of the protective role of some MMR. For 
instance, the 2:1 MMR associated with the SAR consent 
together stable configurations even for enormously high ec- 
centricities, ~ 0.95-0.98 (Gozdziewski et al. 2003 and refer- 
ences therein). This could explain a very large eccentricity 
for planet d and still preserve the system stability by keeping 
the planets away from close encounters. Actually, a modi- 
fication of the relative phase of the planets strongly affects 
the synthetic RV curve and a stable resonance configura- 
tion can be far from being consistent with the RV observa- 
tions. We find that manipulating the values of oJd and of the 
mean anomaly, Md, to get stable configurations is highly 



unfavourable by the RV data (i.e. poor fits are obtained). 
Therefore, this supports the argument arisen in the previ- 
ous paragraph about performing an exploration focused on 
the (Pd - ed) space while leaving to the algorithm the task of 
fitting, without constraints of any sort, the other parameters 
and in particular md, ujd and Md- 

To perform Newtonian orbital fits, Systemic offers dif- 
ferent method such as the Runge-Kutta, Hermite 4 th order 
and Gragg-Bulirsch-Stoer integrators. Fitting a Newtonian 
solution takes longer than a Keplerian model but it assures 
short time scale interactions, relevant for planets c and d, 
are considered. 
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The following step is studying the stability of each dis- 
tinct fit over a period of time related to the time-scale of 
unstable behaviours. For these long-term evolution tests, di- 
rect N-body integrations are applied to the orbital solutions 
considered. We integrate the orbits for at least 1 Myrs us- 
ing the Wisdom-Holman Mapping integrator availa ble in the 
SWIFT software package (jLevison fc Duncanlll994l ). We use 
a time step approximately equal to a twentieth (~ 1/20) of 
the orbital period of the innermost planet. To study close en- 
counters, we use the available Bulirsch-Stoer integrator with 
a tolerance parameter of 10 -9 . We identify each configura- 
tion as being a stable system if orbits stay well bounded over 
an arbitrarily long period of time. 

The results of our analysis are outlined in Figs. [4] and 
[5] Here we label as stable the solutions that survived at 
least for 1 Myrs. Once again, we fix (Pd - &d) and then look 
for the best fit, starting the L-M scheme with initial points 
derived from previously studied configurations. The L-M al- 
gorithm and Amoeba offer a clear representation of the pa- 
rameter space. The dynamical analysis reveals a narrow and 
long band around 3.3 AU and a small island around 3.2 AU 
where good fitness is achieved and stability requirements are 
met. We find a configuration, that we label as stable best- 
fit, which survives for at least 250 Myrs (see Sect. [3] for an 
in-depth examination). Other models scored a better Xred 
but did not preserve stability for the same amount of time. 
Therefore, the stable best-fit seems to lie on the border of 
a chaotic and unstable zone where small changes on the pa- 
rameters of the outermost planet may push the system into 
a strongly chaotic state leading, in some scenarios, even to 
its disruption. 

Figure [4] illustrates how stable configurations do not 
exist in the near neighbourhood of the statistically best- fit; 
smaller quantities for are needed in order for the models 
to retain stability and that increases the value of Xred- 

The top panel of Figure [S] shows the best fits obtained 
during our investigation in terms of the mass for planet d, 
md, and the semi-major axis a<j. The picture makes clear 
how to explain a certain RV amplitude Kd, a bigger mass 
md is required as long as ad increases. 

The bottom panel of Figure [5] illustrates the results of 
our analysis in the semi-major axis-eccentricity plane (ad 
- ed)- The parameters represented are the osculating ele- 
ments at the epoch of the first observation. We show the 
collision line which is defined in terms of semi-major axes 
and eccentricities as a c (l + e c ) = a<i(l + e^). This line de- 
notes the region where the mutual interactions of relatively 
massive companions can rapidly destabilize the configura- 
tion and is calculated for e c = 0.269 and a c — 1.773 AU 
(the values are from our stable best-fit solution, see Table 
[T] later on). Note how the statistically best-fit is positioned 
well over the collision line. We also identify the most relevant 
MMRs between planets c and d, such as the 2:1, 11:5, 9:4, 
7:3, 12:5, 5:2, 8:3, 11:4, 3:1, 10:3 and 7:2. The position of the 
indicated locations has been calculated with respect to the 
values of the stable best-fit. Planets in some resonant config- 
urations, even if under the collision line, exchange angular 
momentum rapidly; their eccentricities are quickly pumped 
and that may lead again to instabilities and self-disruptions. 
In particular, we have found models that show a stable and 
bounded evolution for many Myrs before the unstable be- 
haviours are manifested. On the contrary, other resonant 
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Figure 4. The outcome of the simulations for the HD 181433 
planetary system in terms of the statistical goodness of the fits 
and the orbital period of the outermost planet. The statistically 
best-fit is indicated with an inverted triangle (blue), the stable 
best-fit is denoted with a circle (red). With triangles (green) we 
represent unstable configurations while with squares (black) we 
refer to models stable for at least 1 Myrs. Stable configurations 
do not exist in the near neighbourhood of the statistically best- 
fit, the deep minimum with stable models represents the region 
where a 5:2 MMR configuration is possible. 

configurations, such as the 5:2 and 7:2, are observed to re- 
tain stability even for values over the collision line. 

Table [1] reports the determined set of orbital elements 
for the stable best-fit. For each planet, we list period (P), 
time of periastron passage (T peri ), eccentricity (e), argu- 
ment of pericentre (cj), semi-amplitude (K), minimum mass 
(msini) and semi-major axis (a); we indicate also the stellar 
offset ( V). This model has xled = 4.96, an rms scatter of 
1.36 m/s and the expected jitter of HD 181433 is 1.19 m/s. 
Figure|6]displays the RV data fitted to this model along with 
the residuals. The right panel of figure [3] shows the orbital 
configuration of the system, this time the orbits of planets 
c and d do not cross each other. 

Since the stable best-fit is found in an active region, 
rather than estimating an uncertainty on each parameter, 
we think Figs. [4] and [5] are more useful in visualizing the 
results of the dynamical study and highlight what is plau- 
sible to expect from new observations. If we compare our 
results with what has already been published for this plan- 
etary systenjf], we find that the parameters of planet b and 
c are confirmed to be already well constrained with just K c 
not compatible within the 3 a. For planet d, all the elements 
are found within the 3 a from the original conclusion. How- 
ever, it is worth to underline how to explain the very large 
eccentricity of the third planet and to retain a good fit to 
the present data, the uncertainity on the location of planet 
d reduces drammatically to the narrow band where the 5:2 
MMR is possible. Hence, this supports how a dynamical 
study can be fundamental in interpreting observations, pro- 

4 iBouchv et al.l (2009) do not report directly the uncertainties for 
the masses and semi-major axes, in this case we considerer what 
is available on http://exoplanets.org 
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Table 1. Orbital and physical parameters of the stable best- fit found for the HD 181433 planetary system. 
The osculating elements are given for the epoch of the first observation BJD 2452797.8654. 



Parameter 


HD 181433 b 


HD 181433 c 


HD 181433 d 


P (days) 


9.37459 


975.41 


2468.46 


T peri (BJD-2450000) 


2788.9185 


2255.6235 


1844.4714 


e 


0.38840 


0.26912 


0.46626 


w (°) 


202.039 


22.221 


319.129 


K (m/s) 


2.57 


14.63 


9.41 


msini (mjup) 


0.02335 


0.65282 


0.52514 


msini (mjj) 


7.4 


207.5 


166.9 


a (AU) 


0.08013 


1.77310 


3.29347 


V (m/s) 

rms (m/s) 

Y 2 , 
*-red 




40212.846 
1.36 
4.96 





> 





Julian Days 



2452800 2453150 2453500 2453850 
Time (BJD) 



Figure 6. The Newtonian 3-planet stable best-fit model and residuals periodogram for the HD 181433 RV data. 



during a self-consistent model compatible with the data and 
giving substantial constraints on the orbital parameters. 

The data do not offer any possibility of constraining the 
orbital inclinations. The Newtonian model cannot be partic- 
ularly improved because, aside from the fact the signal of the 
outer planet is not well sampled, we need to wait for secular 
timescales before the variations in i can be spotted via the 
RV method. In fact, we notice that for planets GJ 876 b 
and c which have the strongest mutual gravitational inter- 
actions, more than 11 years of observations (corresponding 
to more than 60 orbits of the outer planet) were used to 
give a reasonable estimate of the inclinations (|Correia et al.l 
l2010t) . 



2.3 Additional planets? 

Finally, following the argument that the proximity of the 
best fit to th e collision line may indica te the presence of fur- 
ther planets (jGozdziewski et al.ll2008l ). we aim to search for 
4-planet Newtonian solutions. The periodogram of the resid- 
uals to the 3-planet solution, in Fig. displays no strong 
peaks that would support the evidence for additional plan- 
ets in the system. Apart from more distant companions, 



in the inner region of the system a terrestrial planet can 
only survive if located between planets b and c. In fact, 
already planet b is found in the proximity of the parent 
star and the area between planets c and d is dominated 
by the strong interactions that interest the two giant plan- 
ets in eccentric orbits. The existance of this last planet 
would support the "packe d planetary systems" hypothesis 
jBarnes fc Ravmondlfiooi ). 



The present data do not allow making any supposition 
about possible outer planets (for example an object at 7 AU 
would have an orbital period of around 7700 days). On the 
other hand, with a super-Earth in the stable zone between 
planets b and c the fit improves. However, this signal would 
be at the noise level with an F-test of the order of 30%. 
The F-test indicates the probability that a planetary model 
would produce a signal similar to the one due just to noise 
fluctuations in the data (Marcy et al. 2005 and references 
therein), so additional observations are required to investi- 
gate on the presence of a super-Earth or less massive planet 
in this stable region. 
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Figure 5. The best fits obtained for the HD 181433 planetary 
system. Top panel: In terms of the mass and semi-major axis of 
the outermost planet. Bottom panel: In terms of the eccentricity 
and semi- major axis of planet d; the collision line is depicted, 
the nominal positions of the most relevant MMR are also labelled 
and marked by dashed lines. The statistically best-fit is indicated 
with an inverted triangle (blue), the stable best-fit is denoted 
with a circle (red). With triangles (green) we represent unstable 
configurations while with squares (black) we refer to models stable 
for at least 1 Myrs. The size of each symbol is proportional to its 
Xred' ' ,e ' smauer symbols indicate better fits. The 5:2 and 7:2 
MMRs retain stability even for values of over the collision line. 



3 LONG-TERM BEHAVIOUR OF THE 
STABLE BEST-FITTING 
CONFIGURATIONS 

Because of the proximity of the two outermost planets, the 
system cannot be stable unless a resonant mechanism is 
present to avoid close encounters. In this Section we aim 
to deepen the study of the stable best-fit configuration as 
well as investigating the evolution of the orbital elements, 
the secular resonant arguments and critical angles of some 
particular configurations consistent with the RV observa- 
tions. 

For the stable best-fit, Fig. [8] shows in the subsequent 
panels the time evolution of the semi-major axes and of the 
eccentricities. Moreover, a secular resonant angle and a crit- 
ical argument of the 5:2 MMR are also illustrated. Specifi- 



Figure 7. Pcriodogram of the residuals to the 3-planct solution 
for HD 181433. It does not display any strong peak that would 
support the evidence for additional planets in the system. The 
three horizontal lines represent, from top to bottom, the 0.1 %, 
1.0 % and 10.0 % FAP levels, respectively. 



cally, the top-left panel in Figure[S]highlights that for the 250 
Myrs of the numerical integration the apocentre of planet c 
and the pericentre of planet d share the same region. If we 
get a close-in view of the situation, we notice that actually 
they never cross each other. In particular, the pericentre of 
d is internal to the apocentre of c. The former approaches 
the value of ad around every 50000 years. It is probably a 
resonance that, protecting the companions from close en- 
counters, allows the stability of the system. The top-right 
panel illustrates the relative large range in which e c and 
ed evolve. The peak-to-peak amplitude is covered in around 
2500 years only. e c moves in the interval 0.17-0.52 while ed 
in the range 0.17-0.50. The present eccentricities fall in the 
middle of these intervals indicating that the system has been 
snapped in a statistically quite probable state. Also, such 
a large range reminds that in multiple-planet systems the 
orbital eccentricities can vary considerably through secular 
interactions on timescales that are long compared to ob- 
servational baselines but short compared to the age of the 
systems. Therefore, when doing statistical studies on exo- 
planetary systems the planetary orbits should normally be 
described by a complete distribution of values for the eccen- 
tricities rather than just by the present quantities (see also 
Adams & Laughlin 2006). This model is not observed to be 
in SAR. The bottom-left panel indicates the time evolution 
of the secular argument uj c + 0Jd which alternates librations 
with circulations. Besides, we find that 5rid—2n c « — 3.4° /yr 
indicating the proximity of the 5c:2d MMR, therefore here 
we have a scenario similar to the Jupiter-Saturn case in the 
Solar System. We have studied the resonant arguments of 
this MMR and found that the angle 5Xd — 2A C — 3uud librates 
around 180° with a semiaplitude of about 110°. Thus, this 
configuration is seen to be locked in a MMR, the bottom- 
right panel illustrates how this critical angle evolves with 
time. 

The mass of planet b is negligible with respect c and 
d so we can assume the dynamics of the two giants is not 
disturbed much by the presence of the rocky planet close by. 
Then, we study some possible configurations with planets c 
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Figure 8. Evolution of some orbital elements for the stable best-fit configuration. Top panels: The time evolution of the semi-major axes 
and of the eccentricities are depicted, the solid line (black) indicates planet b, the dashed line (red) denotes planet c and the dotted line 
(blue) represents planet d. e c moves in the interval 0.17-0.52 while in the range 0.17-0.50. Bottom panels: A secular resonant angle 
and a critical argument of the 5:2 MMR are illustrate. 



and d near MMR. We do not find any plausible (Xred ^ 6.03) 
stable solution that would correspond to the 2:1, 11:5, 9:4, 
7:3, 8:3, 11:4 and 10:3 MMRs. In particular, the configura- 
tions nMMR 11:5, 9:4 and 7:3 seem preferred by the data 
but the excessive pumping of the eccentricities causes close 
encounters and planetary scatterings which do not favour 
stability. On the other hand, outer MMRs are possible be- 
cause here in particular the planets are more spread and 
collisions can be avoided. 

We compute the evolution of the orbital elements for 
fits corresponding to MMRs 12:5, 7:3, 3:1 and 7:2 which our 
simulations have demonstrated to preserve stability for at 
least 40 Myrs. The results are illustrated in Figs. lit flTJl fTTI 
and [12] and show the complexity of the possible dynamical 
behaviours of the HD 181433 system that are consistent with 
the RV observations. 

The case we show nMMR 12:5 has Xred = 5-91 an d rms 
scatter of 1.46 m/s with 12rid — 5n c ~ 1.8° /yr. This model is 
in SAR with lj c — uud librating around 0° with a semiampli- 
tude of about 45°. Moreover, we have found some critical 
arguments of the MMR to alternate librations with circula- 



tions implying the resonance excites a chaotic configuration. 
The time evolution of the critical angles 12Ad — 5A C — 7u c 
and 12Ad — 5A C — 5cj d — 2u) c are illustrated in Fig. [5] 

The scenario nMMR 7:3 has Xred = 5-89 and rms scat- 
ter of 1.46 m/s with 7na — 3n c ~ 8.7° /yr. This model is ob- 
served to be in SAR with the critical angle lj c — uid librating 
around 0° with a semiamplitude of about 40° meaning that 
their periastrons are aligned. We have studied the resonant 
arguments of the 7:3 MMR and found that three angles i.e. 
7Ad — 3A C — Auid, 7Ad — 3A C — 3oJd— ^ c and 7Ad — 3A C — 2u>d — 2uj c , 
alternate librations with circulations. This indicates that the 
configuration is close to the resonance separatrices. Fig. [10] 
illustrates how the critical angle 7Ad — 3A C — 3ojd — oj c evolves 
with time. 

The case near the low-order MMR 3:1 has Xred = 6-03 
and rms scatter of 1.47 m/s with 3rid — n c ~ 0.6° /yr. This 
model is in SAR with ui c — u>d librating around 180° with 
a semiamplitude of about 110° implying that this time the 
periastrons of planets c and d are anti-aligned. The data 
points that diverge from the periodic signal (Fig. [11] top- 
right panel) represent the instants when ed gets close to 
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Figure 9. Evolution of some orbital elements for the nMMR 12:5 configuration. Top-left: The solid line (black) indicates planet b, the 
dashed line (red) denotes planet c and the dotted line (blue) represents planet d. e c moves in the interval 0.14-0.31 while in the range 
0.11-0.31. This model is in SAR. 
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Figure 10. Evolution of some orbital elements for the model nMMR 7:3. Left Panel: The solid line (black) indicates planet b, the 
dashed line (red) denotes planet c and the dotted line (blue) represents planet d. e c moves in the interval 0.17-0.30 while in the range 
0.15-0.30. This model is observed to be in SAR. 



be null and so the argument of pericentre, uid, is not well 
defined. Furthermore, for the critical arguments of the MMR 
we find that librations alternate with circulations indicating 
a chaotic zone spanned by overlapping resonances. The time 



evolution of the critical angles 3A<j — A c — U)d — uj c and 3Ad — 
A c — 2uj c are also illustrated in Fig. 1111 

The scenario nMMR 7:2, illustrated in Fig. [12] has 
Xred = 5-65 and rms scatter of 1.45 m/s with 7n<j — 2n c w 
— 1.7° /yr. This model is not seen to be in SAR. However, 
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Figure 11. Evolution of some orbital elements for the model nMMR 3:1. Top-left: The solid line (black) indicates planet b, the dashed 
line (red) denotes planet c and the dotted line (blue) represents planet d. e c moves in the interval 0.23-0.37 while in the range 0.0-0.27. 
This model is in SAR. 



we find it to be locked in the MMR as the critical argument 
7Ad — 2A C — 5uJd librates around 180° with a semiamplitude 
of about 85°. It is worth notice how this mechanism is ca- 
pable of pumping e c from 0.09 to 0.47 in less than 5000 
years. This configuration is located over the collision line, it 
is the resonance that prevents close encounters and provides 
long-term orbital stability to the system. 

Considering all the the configurations studied, the be- 
haviour of planet b seems to be unrelated to the two giant 
companions as the amplitude of the eccentricity signal of b 
appears to be unaffected even in the cases in which c and 
d are trapped in a resonance. Moreover, it seems unrealistic 
that planet b can be involved in a p:q:r MMR with the outer 
two planets: b is too far away from them and, in addition, it 
is mainly influenced by General Relativistic and tidal effects 
(a discussion of these mechanisms goes beyond the aim of 
this paper). 

The synthetic RV curves for the Keplerian best- fit, sta- 
ble best-fit and the models nMMRs 12:5, 7:3 and 7:2 are 
illustrated in Fig. 1131 The period through the year 2017 is 
covered. It is difficult to distinguish the curves from each 
other in the time range covered by the observations, but 



at the time of writing the difference can be spotted. How- 
ever, the signal of the Keplerian best-fit will diverge evi- 
dently from the stable best-fit only in February 2013 (~ JD 
2,456,340). 



4 DISCUSSION AND CONCLUSION 

In this paper, our efforts have been focused on finding a 
plausible solution to the available RV data for the planetary 
system of HD 181433. In our investigation we have included 
an analysis of the long-term evolution of the system and 
the results support the thesis the dynamics is an important 
observable that has to be taken into account with the same 
priority as the RV observations. 

The story with HD 181433 does not differ from the one 
of many other multi-planet systems found on the edge of 
long-term dynamical stability. The dynamical modelling of 
the RV with stability constraints offers precious information 
on the dynamical architecture of the putative planetary con- 
figurations. The stability criterion becomes a fundamental 
tool which provides limits to the physical and orbital ele- 
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Figure 12. Evolution of some orbital elements for the model nMMR 7:2. Left Panel: The solid line (black) indicates planet &, the 
dashed line (red) denotes planet c and the dotted line (blue) represents planet d. e c moves in the interval 0.09-0.47 while ed in the range 
0.32-0.53. The critical argument illustrated in the right panel librates around 180° with a semiamplitude of about 85°. 
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Figure 13. Synthetic RV curves for HD 181433. In every panel, the straight line (black) indicates the Keplerian best-fit while the dotted 
line (red) represents (clockwise starting from the top-left) the stable best-fit, the configuration studied near the MMR 12:5, the model 
nMMR 7:2 and the scenario nMMR 7:3. Data points (blue) are plotted with error bars (and V ff se t is the one calculated for the stable 
best-fit). The dashed (blue) vertical line indicates when the last data point was taken. The straight (black) vertical line at JD 2,455,760 
corresponds to the calendar date 2011 July 17 when the difference between most of the curves can be identified. 
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ments of the planets when data do not cover completely the 
longest orbital period. 

Our investigation leads to a Newtonian model for HD 
181433 stable for at least 250 Myrs. The sol u tion i s compat- 
ible with what was found by iBouchv et al l (2009), but our 
analysis strongly diminish the uncertainity on the location of 
planet d to the exiguous band where the 5:2 MMR is possi- 
ble and stability is preserved. This seems the only plausible 
way to explain a very large eccentricity for the outermost 
planet, a quality which must be met in order to hold a good 
fit to the present data. In general, we can say that when an 
unstable high eccentric solution is found for a multiplanet 
system, the study of resonances may lead to the finding of 
a reasonable stable solution. By doing so, it is possible to 
constrain with high confidence the orbital period of the out- 
ermost (poorly sampled) planet well before sufficent data, 
covering several orbital periods, become available. 

Apart the 5:2 MMR, the orbital evolution of the two gi- 
ant companions is confined to a zone spanned by a number of 
other low-order two-body MMRs. We have studied different 
plausible stable configurations for the planetary system and, 
in particular, we have illustrated the behaviours caused by 
secular apsidal resonances and mean motion resonances. We 
have also found that at the time of writing with new data 
points it is definitely feasible to refine the circle of likely 
scenarios. 

Furthermore, given the strong gravitational interactions 
between the two giant planets, a self-consistent N-body 
model for the RV data will help in estimating the incli- 
nation of the planetary orbits and of the physical masses 
(the RV method returns just the minimum masses for the 
planets). If it is not possible to use in synergy other meth- 
ods e.g. astrometry, transits, it may be necessary to observe 
around 50 full orbits of planet d i.e. 300 years of RV data, 
to strongly constrain the orbital inclinations. As the large 
values already observed for the eccentricities may have been 
trigged by mechanisms which influence also the orbital incli- 
nations (e.g. Libert & Tsiganis 2009), a considerable value 
for their mutual inclination may be expected. 

We can calculate the habitable zone (HZ) orbital dis- 
tance dhat, to be defined as the distance where a planet 
would receive the same insolation as the Earth: aw a b = \ I -r 2 - 

V L <3 

AU. For HD 181433, we find a hab = 0.55 AU (P hab ~ 170 
days) which is in the region between planets b and c. Our 
simulations, run for 40 Myrs, shows an Earth-size planet in 
the HZ (and in eccentric orbit) can retain stability indeed. 
The existence of a planet in this zone, not only fills an empty 
gap in the system, but would also offer a harbour for life. 
Therefore, this hypothetical planet becomes interesting for 
a double reason. Additional observations are required to in- 
vestigate on the presence of further bodies and in particular 
on the existence of a terrestrial planet in the habitable zone 
of HD 181433. 

The planetary system of HD 181433 offers a wide range 
of challenges that can go from understanding its real physi- 
cal architecture to the study of potentially habitable worlds. 
Further observations can confirm the results illustrated in 
this paper, improve our understanding of the dynamical 
structure of this system and, in general, give additional in- 
sights into the study of the dynamics of planetary systems. 
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